Port Everglades Coral Survey Data Analysis
R. Cunning
2025-04-21
# Load packages
library(sf)
library(xml2)
library(tidyverse)# Define genus level taxon groups (plus one family FAVI)
taxon_groups <- list(
PORI = c("PPOR", "PFUR", "PDIV", "PAST", "PORI"),
ORBI = c("OFAV", "OANN", "OFRA", "ORBI"),
FAVI = c("CNAT", "DLAB", "PSTR", "PCLI", "MARE", "FAVI"),
AGAR = c("AFRA", "AAGA", "AHUM", "ALAM", "AGAR"),
MADR = c("MAUR", "MSEN", "MDEC", "MPHA", "MADR"),
SOLE = c("SHYA", "SBOU", "SOLE"),
SCOL = c("SLAC", "SCUB", "SCOL"),
SIDE = c("SSID", "SRAD", "SIDE"),
MYCE = c("MFER", "MLAM", "MALI", "MYCE")
)
# Convert to lookup tibble
taxon_lookup <- enframe(taxon_groups, name = "taxon_group", value = "taxon") %>%
unnest(taxon)
# Define juvenile family level taxon groups (following DRM survey convention)
taxon_groups_juv <- list(
MUSS = c("ISIN", "ISOP", "MANG", "MYCE", "SCOL", "MUSS"),
FAVI = c("FAVI", "FFRA", "MARE"),
MEAN = c("MMEA", "MEAN", "DCYL", "DSTO", "EFAS")
)
# Convert to lookup tibble
taxon_lookup_juv <- enframe(taxon_groups_juv, name = "taxon_group", value = "taxon") %>%
unnest(taxon)1. Import data
1.1. 2024 Shedd Survey
1.1.1. Shedd site metadata
# Site metadata
drm_sitemd <- readxl::read_xlsx("data/2024_shedd_drm/site_metadata.xlsx") %>%
janitor::clean_names() %>%
mutate(site = as.character(drm_site_id)) %>%
select(site, latitude = lat, longitude = lon) %>%
mutate(depth = NA)1.1.2. Shedd coral data
# Adult coral data from main DRM surveys
#adults0 <- read_csv("data/2024_shedd_drm/DRM_broward_corals.csv")
#adults0 %>% filter(team == "Shedd Aquarium")
# Most sites were included in main DRM database for 2024 -- Import these
adultst1t2 <- readxl::read_xlsx("data/2024_shedd_drm/2024ANU_RawCoralDataTransect1and2_Shedd.xlsx") %>%
janitor::clean_names() %>%
filter(subregion == "Broward-Miami", team == "Shedd Aquarium") %>%
select(site, transect_num, species, width, height)
adultst3t4 <- readxl::read_xlsx("data/2024_shedd_drm/2024ANU_RawCoralDataTransect3and4_Shedd.xlsx") %>%
janitor::clean_names() %>%
filter(subregion == "Broward-Miami", team == "Shedd Aquarium") %>%
select(site, transect_num, species, width, height)
# 9 of our PEV sites were removed from DRM database to avoid oversaturating the ares -- Import these separately
removedt1t2 <- readxl::read_xlsx(
"data/2024_shedd_drm/2024_DRM_Broward_RemovedSites_T1-T4_Shedd.xlsx", sheet = "Removed Sites T1-T2") %>%
janitor::clean_names() %>%
select(site, transect_num, species, width, height)
removedt3t4 <- readxl::read_xlsx(
"data/2024_shedd_drm/2024_DRM_Broward_RemovedSites_T1-T4_Shedd.xlsx", sheet = "Removed Sites T3-T4") %>%
janitor::clean_names() %>%
select(site, transect_num, species, width, height)
# Combine all adult coral data for Shedd DRM surveys at PEV
adults0 <- bind_rows(adultst1t2, adultst3t4, removedt1t2, removedt3t4)
# Convert adult data to long format
adults_long <- adults0 %>%
mutate(max_width_cm = pmax(width, height, na.rm = TRUE)) %>%
mutate(max_width_cm = as.character(max_width_cm)) %>%
mutate(site = str_remove(site, "^AA")) %>%
select(site, transect_num, taxon = species, max_width_cm) %>%
drop_na(taxon) # DROPS when taxon is blank, this is when no corals >4cm were observed
# Juveniles from main DRM surveys (family-level tallies)
# drmjuv <- read_csv("data/2024_shedd_drm/DRM_broward_juveniles.csv") %>%
# filter(team == "Shedd Aquarium") %>%
# mutate(site = str_remove(site, "^AA")) %>%
# select(team, site, transect_num, ends_with("ct")) %>%
# rename(MCAV = montastraea_ct, MUSS = mussinae_ct, FAVI = faviinae_ct, MEAN = meandrinidae_ct)
# Import juvenile counts from main DRM dataset
juv <- readxl::read_xlsx("data/2024_shedd_drm/2024ANU_JuvenileCoralData_Shedd.xlsx") %>%
janitor::clean_names() %>%
filter(subregion == "Broward-Miami", team == "Shedd Aquarium")
# Import juvenile counts from sites that were removed from main DRM dataset
removed_juv <- readxl::read_xlsx("data/2024_shedd_drm/Shedd_removed_sites_Juveniles_2024.xlsx") %>%
janitor::clean_names() %>%
# Missing values in count data should be zero counts (unique to this datasheet from FWC)
mutate(across(ends_with("_ct"), ~replace_na(., 0)))
# Combine juvenile data
juv0 <- bind_rows(juv, removed_juv) %>%
mutate(site = str_remove(site, "^AA")) %>%
select(site, transect_num, ends_with("ct")) %>%
rename(MCAV = montastraea_ct, MUSS = mussinae_ct, FAVI = faviinae_ct, MEAN = meandrinidae_ct)
# Convert juvenile data to long format
juv_long <- juv0 %>%
pivot_longer(c(MUSS, FAVI, MEAN, MCAV), names_to = "taxon", values_to = "n") %>%
mutate(max_width_cm = "<4") %>%
uncount(n)
# Other juvenile taxa counts from Transects 1 and 2 (DRM 'bonus data')
t1t2bonus <- read_csv("data/2024_shedd_drm/T1_T2_bonus_data.csv") %>%
janitor::clean_names() %>%
mutate(site = replace_na(site, "NA")) %>% # Because one site is called "NA"
mutate(transect_num = parse_number(transect))
t1t2juv <- t1t2bonus %>%
select(site, transect_num, starts_with("small")) %>%
rename_with(~ toupper(gsub("^small_", "", .x)), starts_with("small_"))
t1t2juv_long <- t1t2juv %>%
pivot_longer(3:10, names_to = "taxon", values_to = "n") %>%
mutate(max_width_cm = "<4") %>%
uncount(n)
# Replace site names in t1t2 bonus data with the correct DRM site ID
penipsites <- readxl::read_xlsx("data/2024_shedd_drm/site_metadata.xlsx") %>%
janitor::clean_names()
t1t2juv_long_updated <- t1t2juv_long %>%
left_join(penipsites %>% select(site, drm_site_id), by = "site") %>%
mutate(site = as.character(drm_site_id)) %>%
select(-drm_site_id)
# Combine all data
drm_long <- bind_rows(adults_long, juv_long, t1t2juv_long_updated) %>%
mutate(team = "Shedd Aquarium")
# Check species names
sort(unique(drm_long$taxon))## [1] "AAGA" "ACER" "AFRA" "CNAT" "DLAB" "DSTO" "EFAS" "FAVI" "FFRA" "MALI"
## [11] "MAUR" "MCAV" "MDEC" "MEAN" "MMEA" "MUSS" "OANN" "OFAV" "OFRA" "PAST"
## [21] "PCLI" "PFUR" "PPOR" "PSTR" "SBOU" "SCUB" "SINT" "SRAD" "SSID"write_csv(drm_long, file = "data/processed/drm_2024_long.csv")
# COUNT based on rules
# ✅ Updated Rules Summary for counting from DRM/Shedd data:
# Juvenile taxa (searched for in <4cm size class only):
# "MEAN", "MUSS", "FAVI"
# → these should only ever appear in <4cm, never >4cm, and should not be zero-filled for adults.
# Other juvenile-capable taxa:
# "MCAV", "SSID", "SRAD", "PAST", "PPOR", "SINT", "SBOU", "AAGA", "MAUR"
# → these can be counted in both >4cm and <4cm, but only in <4cm if juveniles were searched on that transect and team.
# Transect-based search rules still apply:
# Transects 1 & 2: all adult taxa always searched. Juvenile search depends on team:
# "Shedd Aquarium" → all juvenile taxa above searched
# others → only MEAN, MUSS, FAVI, MCAV
# Transects 3 & 4:
# only subset of adult taxa searched (adult_taxa_t3t4)
# only juveniles: MEAN, MUSS, FAVI, MCAV
# Step 1: Define size classes
drm_classed <- drm_long %>%
mutate(class = case_when(as.numeric(max_width_cm) >= 4 ~ ">4cm",
max_width_cm == "<4" ~ "<4cm"))
# Step 2: Define species sets
all_taxa <- unique(drm_classed$taxon)
adult_taxa_t3t4 <- c("CNAT", "DSTO", "DLAB", "MMEA", "MANG", "MALI",
"MFER", "MLAM", "PCLI", "PSTR")
juv_only_taxa <- c("MEAN", "MUSS", "FAVI")
juv_both_taxa <- c("MCAV", "SSID", "SRAD", "PAST", "PPOR", "SINT", "SBOU", "AAGA", "MAUR")
all_juv_taxa <- c(juv_only_taxa, juv_both_taxa)
# Step 3: Build search grid per site × transect × team
search_grid <- drm_classed %>%
distinct(site, transect_num, team) %>% # if multiple teams in data, remove value for team
mutate(
searched_taxa_class = pmap(list(transect_num, team), function(transect, team) {
# Helper: define juv taxa allowed for this transect/team
juv_taxa <- if (transect %in% c(1, 2)) {
if (team == "Shedd Aquarium") { # Shedd searched for other juv taxa on T1 and T2, other DRM survey teams did not
all_juv_taxa
} else {
c(juv_only_taxa, "MCAV")
}
} else {
c(juv_only_taxa, "MCAV")
}
# Adults always searched in 1 & 2, subset in 3 & 4
adult_taxa <- if (transect %in% c(1, 2)) {
setdiff(all_taxa, juv_only_taxa) # exclude juv-only taxa
} else {
adult_taxa_t3t4
}
# Build grid
bind_rows(
expand_grid(taxon = adult_taxa, class = ">4cm"),
expand_grid(taxon = juv_taxa, class = "<4cm")
)
})
) %>%
unnest(searched_taxa_class)
# Step 4: Count observations
counts <- drm_classed %>%
group_by(site, transect_num, team, taxon, class) %>%
summarize(n = n(), .groups = "drop")
# Step 5: Join with grid and fill in zeros where appropriate
final_counts <- search_grid %>%
left_join(counts, by = c("site", "transect_num", "team", "taxon", "class")) %>%
mutate(n = replace_na(n, 0))
write_csv(final_counts, file = "data/processed/drm_2024_counts.csv")
# AGGREGATE COUNT DATA
# Aggregate taxa
final_counts_ag <- final_counts %>%
left_join(taxon_lookup, by = "taxon") %>%
mutate(taxon = coalesce(taxon_group, taxon)) %>%
select(-taxon_group) %>%
group_by(across(-n)) %>%
summarize(n = sum(n), .groups = "drop")
write_csv(final_counts_ag, file = "data/processed/drm_2024_counts_ag.csv")1.2. 2017 Dial Cordy Recon Survey
1.2.1. DCA site metadata
# Site metadata
# Read in site coordinates
dca_sitemd0 <- readxl::read_xlsx("data/2017_dca/Recon_Site_Coordinates_Extracted.xlsx") %>%
janitor::clean_names()
# All sites have start and end coordinates...
# Tidy and Calculate midpoint per transect
dca_sitemd <- dca_sitemd0 %>%
mutate(
depth = abs(as.numeric(depth)),
across(c(latitude, longitude), as.numeric)
) %>%
group_by(site = transect) %>%
summarize(
latitude = mean(latitude, na.rm = TRUE),
longitude = mean(longitude, na.rm = TRUE),
depth = mean(depth, na.rm = TRUE),
.groups = "drop"
)1.2.2. DCA coral data
# Read in survey data
dca0 <- readxl::read_xlsx("data/2017_dca/Compiled_DCA_RECON_Belt_data.xlsx") %>%
janitor::clean_names()
dca <- dca0 %>%
select(1:18) %>%
rename(site = site_name) %>%
mutate(site = factor(site)) %>%
select(site, taxon = coral_species, max_width_cm = max_size_cm)
# Adjust/corrects species IDs
sort(unique(dca$taxon))## [1] "AAGA" "ACER" "AFRA"
## [4] "AGA SP" "AHUM" "ALAM"
## [7] "CNAT" "CORAL" "Cup Coral"
## [10] "DLAB" "DSTO" "EFAS"
## [13] "LCUC" "MAD SP" "MADSP"
## [16] "MALI" "MARE" "MCAV"
## [19] "MDEC" "MLAM" "MMEA"
## [22] "MPHA" "Mycetophyllia spp." "MYCSP"
## [25] "OANN" "ODIF" "OFAV"
## [28] "OFAV\\" "PAST" "PCLI"
## [31] "PFUR" "PPOR" "PSTR"
## [34] "SBOU" "Scolymia Spp" "SCUB"
## [37] "Sid SP" "SID SP" "SID SP."
## [40] "SIDSP" "SINT" "SRAD"
## [43] "SSID"dca <- dca %>%
mutate(taxon = case_when(
taxon == "AGA SP" ~ "AGAR",
taxon == "LCUC" ~ "HCUC",
taxon %in% c("MYCSP", "Mycetophyllia spp.") ~ "MYCE",
taxon == "OFAV\\" ~ "OFAV",
taxon %in% c("MAD SP", "MADSP") ~ "MADR",
taxon == "Scolymia Spp" ~ "SCOL",
taxon %in% c("SIDSP", "Sid SP", "SID SP.", "SID SP") ~ "SIDE",
TRUE ~ taxon
))
sort(unique(dca$taxon))## [1] "AAGA" "ACER" "AFRA" "AGAR" "AHUM" "ALAM"
## [7] "CNAT" "CORAL" "Cup Coral" "DLAB" "DSTO" "EFAS"
## [13] "HCUC" "MADR" "MALI" "MARE" "MCAV" "MDEC"
## [19] "MLAM" "MMEA" "MPHA" "MYCE" "OANN" "ODIF"
## [25] "OFAV" "PAST" "PCLI" "PFUR" "PPOR" "PSTR"
## [31] "SBOU" "SCOL" "SCUB" "SIDE" "SINT" "SRAD"
## [37] "SSID"# Filter out unidentified corals
dca <- dca %>%
filter(!taxon %in% c("CORAL", "Cup Coral"))
# Write long data to file
write_csv(dca, file = "data/processed/dca_2017_long.csv")
# Convert to count data
# Add explicit zeros for any taxon/size class missing at any site
dca_counts <- dca %>%
mutate(class = ifelse(max_width_cm >= 4, ">4cm", "<4cm")) %>%
count(site, taxon, class) %>%
complete(site, taxon, class = c(">4cm", "<4cm"), fill = list(n = 0))# %>%
# # Don't create zeros for MEAN/MUSS/FAVI adults, since these IDs only applied to juv
# filter(!(taxon %in% c("MEAN", "MUSS", "FAVI") & class == ">4cm" & n == 0))
write_csv(dca_counts, file = "data/processed/dca_2017_counts.csv")
# Aggregate count data based on taxonomic groups defined above
dca_counts_ag <- dca_counts %>%
left_join(taxon_lookup, by = "taxon") %>%
mutate(taxon = coalesce(taxon_group, taxon)) %>%
select(-taxon_group) %>%
group_by(across(-n)) %>%
summarize(n = sum(n), .groups = "drop")
# Further aggregate juvenile counts to family (following DRM methods)
dca_counts_ag <- dca_counts_ag %>%
left_join(taxon_lookup_juv, by = "taxon") %>%
mutate(
taxon = if_else(class == "<4cm" & !is.na(taxon_group), taxon_group, taxon)
) %>%
select(-taxon_group) %>%
group_by(across(-n)) %>%
summarize(n = sum(n), .groups = "drop")
write_csv(dca_counts_ag, file = "data/processed/dca_2017_counts_ag.csv")1.3. 2023 Tetra Tech Impact Survey
1.3.1. TT site metadata
# Site metadata
tt_sitemd <- readxl::read_xlsx("data/2023_tt/Impact site tracking.xlsx", skip = 1) %>%
janitor::clean_names()
tt_sitemd <- tt_sitemd %>%
mutate(site = transect_name,
latitude = as.numeric(actual_start_y),
longitude = as.numeric(actual_start_x)) %>%
select(site, latitude, longitude)
# Many sites missing coords in sheet.... whats up with that
tt_sitemd <- drop_na(tt_sitemd, longitude)1.3.2. TT coral data
tt0 <- readxl::read_xlsx("data/2023_tt/Impact Raw Data 05 31 2024.xlsx") %>%
janitor::clean_names()
tt <- tt0 %>%
select(site = transect_name, depth_ft_start,
taxon = id_abbrev, coral_length_cm, coral_width_cm) %>%
filter(taxon != "Xesto") %>%
mutate(taxon = toupper(taxon)) %>%
mutate(across(c(coral_length_cm, coral_width_cm), as.numeric)) %>%
mutate(site = factor(site))
tt <- tt %>%
mutate(max_width_cm = pmax(coral_length_cm, coral_width_cm)) %>%
select(site, taxon, max_width_cm)
# Check taxa names
sort(unique(tt$taxon))## [1] "?" "AAGA" "ACER" "AFRA" "AFRAG" "ALAM"
## [7] "ASP" "ASP." "CNAT" "DLAB" "DSTO" "EFAS"
## [13] "FFRA" "HCUC" "ID-ABBREV" "MCAV" "MCAV?" "MDEC"
## [19] "MHEARD" "MMEA" "MSEN" "MSP." "MUSSID" "MYALI"
## [25] "MYFER" "MYLAM" "NONE" "OFAV" "OFR" "OFRA"
## [31] "OROB" "PAME" "PAST" "PCLI" "PCLI?" "PDIV"
## [37] "PPOR" "PSP" "PSP." "PSTR" "SBOU" "SCUB"
## [43] "SHYA" "SINT" "SLAC" "SRAD" "SSID" "SSP."
## [49] "STOK"# Filter out unidentified taxa
tt <- tt %>%
filter(!taxon %in% c("?", "ID-ABBREV", "NONE", "MHEARD"))
# Adjust/corrects species IDs
tt <- tt %>%
mutate(taxon = case_when(
taxon == "AFRAG" ~ "AFRA",
taxon %in% c("ASP", "ASP.") ~ "AGAR",
taxon == "MCAV?" ~ "MCAV",
taxon == "MSP." ~ "MADR",
taxon == "MUSSID" ~ "MUSS",
taxon == "MYALI" ~ "MALI",
taxon == "MYFER" ~ "MFER",
taxon == "MYLAM" ~ "MLAM",
taxon == "OFR" ~ "OFRA",
taxon == "PCLI?" ~ "PCLI",
taxon %in% c("PSP", "PSP.") ~ "PORI",
taxon == "SSP." ~ "SIDE",
taxon == "STOK" ~ "DSTO",
TRUE ~ taxon
))
sort(unique(dca$taxon))## [1] "AAGA" "ACER" "AFRA" "AGAR" "AHUM" "ALAM" "CNAT" "DLAB" "DSTO" "EFAS"
## [11] "HCUC" "MADR" "MALI" "MARE" "MCAV" "MDEC" "MLAM" "MMEA" "MPHA" "MYCE"
## [21] "OANN" "ODIF" "OFAV" "PAST" "PCLI" "PFUR" "PPOR" "PSTR" "SBOU" "SCOL"
## [31] "SCUB" "SIDE" "SINT" "SRAD" "SSID"# Write long data to file
write_csv(tt, file = "data/processed/tt_2024_long.csv")
# Count
# Add explicit zeros for any taxon/size class missing at any site
tt_counts <- tt %>%
mutate(class = ifelse(max_width_cm >= 4, ">4cm", "<4cm")) %>%
count(site, taxon, class) %>%
complete(site, taxon, class = c(">4cm", "<4cm"), fill = list(n = 0)) %>%
# Don't create zeros for MEAN/MUSS/FAVI adults, since these IDs only applied to juv
filter(!(taxon %in% c("MEAN", "MUSS", "FAVI") & class == ">4cm" & n == 0))
write_csv(tt_counts, file = "data/processed/tt_2024_counts.csv")
# Aggregate count data based on taxonomic groups defined above
tt_counts_ag <- tt_counts %>%
left_join(taxon_lookup, by = "taxon") %>%
mutate(taxon = coalesce(taxon_group, taxon)) %>%
select(-taxon_group) %>%
group_by(across(-n)) %>%
summarize(n = sum(n), .groups = "drop")
# Further aggregate juvenile counts to family (following DRM methods)
tt_counts_ag <- tt_counts_ag %>%
left_join(taxon_lookup_juv, by = "taxon") %>%
mutate(
taxon = if_else(class == "<4cm" & !is.na(taxon_group), taxon_group, taxon)
) %>%
select(-taxon_group) %>%
group_by(across(-n)) %>%
summarize(n = sum(n), .groups = "drop")
write_csv(tt_counts_ag, file = "data/processed/tt_counts_ag.csv")1.4. Habitat classification polygons
# --- STEP 1: Load KML Polygons ---
polygons <- st_read("data/Habitat classifications.kml") # update path as needed## Reading layer `Habitat classifications' from data source
## `/Users/rosscunning/Projects/PENIP/data/Habitat classifications.kml'
## using driver `KML'
## Simple feature collection with 190 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: -80.11618 ymin: 25.97477 xmax: -80.06143 ymax: 26.25943
## z_range: zmin: 0 zmax: 0
## Geodetic CRS: WGS 84# --- STEP 2: Extract Attributes from HTML Description ---
extract_attrs <- function(desc) {
if (is.na(desc) || desc == "") {
return(tibble(Habitat = NA, Type = NA, Modifier = NA, Region = NA, Type2 = NA))
}
html <- read_html(desc)
rows <- xml_find_all(html, "//table//table//tr")
keys <- rows %>% xml_find_all(".//td[1]") %>% xml_text(trim = TRUE)
vals <- rows %>% xml_find_all(".//td[2]") %>% xml_text(trim = TRUE)
n <- min(length(keys), length(vals))
named_vals <- set_names(vals[1:n], keys[1:n])
tibble(
Habitat = named_vals[["Habitat"]],
Type = named_vals[["Type"]],
Modifier = named_vals[["Modifier"]],
Region = named_vals[["Region"]],
Type2 = named_vals[["Type2"]]
)
}
# Apply function and combine with spatial geometries
attrs <- map_dfr(polygons$Description, extract_attrs)
polygons_clean <- bind_cols(polygons %>% select(-Description), attrs)
# --- STEP 3: Prepare Site Coordinate Data ---
# Get all site coordinates
allsitemd <- bind_rows(.id = "dataset",
drm = drm_sitemd,
dca = dca_sitemd,
tt = tt_sitemd
)
points <- st_as_sf(allsitemd, coords = c("longitude", "latitude"), crs = 4326)
# --- STEP 4: Validate Geometry and Match CRS ---
polygons_clean <- polygons_clean %>%
st_zm(drop = TRUE, what = "ZM") %>%
st_make_valid() %>%
st_transform(st_crs(points))
sf_use_s2(FALSE) # prevent s2 geometry issues
# --- STEP 5: Spatial Join ---
joined <- st_join(points, polygons_clean, join = st_within)
joined_df <- joined %>%
mutate(longitude = st_coordinates(.)[,1],
latitude = st_coordinates(.)[,2]) %>%
st_drop_geometry()
allsitemd <- joined_df
# Visualize habitat classifications
(polyplot <- polygons_clean %>%
#filter(Type != "Sand") %>%
ggplot() +
geom_sf(aes(fill = Type), color = "black", size = 0.2, alpha = 0.6) +
scale_fill_brewer(palette = "Set3", na.value = "gray80") +
theme_minimal() +
labs(title = "Habitat Polygons Colored by Type", fill = "Type") +
theme(legend.position = "right") +
xlim(-80.11, -80.079) +
ylim(26.065, 26.11))
# 'Sand' overlaps some of the other polygons instead of just surrounding them...
# If a point is classified as Sand AND something else, remove Sand...
multiclass <- allsitemd %>%
group_by(site) %>%
filter(n() > 1)
# polyplot +
# geom_point(data = multiclass, aes(x = longitude, y = latitude),
# pch = "X", inherit.aes = FALSE)
allsitemd_clean <- allsitemd %>%
group_by(site) %>%
filter(!(Type == "Sand" & n() > 1)) %>%
ungroup()2. Subset comparable survey data
2.1. Find habitat and spatial overlap across surveys
# Plot all sites and select subset based on overlap
polyplot +
geom_point(data = allsitemd_clean,
aes(x = longitude, y = latitude, shape = dataset),
inherit.aes = FALSE, alpha = 0.6) +
#geom_text(data = allsitemd_clean, aes(x = longitude, y = latitude, label = site)) +
scale_shape_manual(values = c(2, 3, 4))
# Select just Nearshore Ridge Complex, Inner Reef, and Middle Reef, because 2024 surveys did not include outer reef
overlap_lat <- allsitemd_clean %>%
filter(latitude > 26.08, latitude < 26.11,
Type %in% c("Nearshore Ridge Complex", "Inner Reef", "Middle Reef", "Artificial"))
overlap_nolat <- allsitemd_clean %>%
filter(Type %in% c("Nearshore Ridge Complex", "Inner Reef", "Middle Reef", "Artificial"))
#anti_join(overlap_nolat, overlap_lat)
# Choose whether to also apply a latitudinal filter since DRM has a few sites further south than DCA/TT
overlap <- overlap_nolat
# Plot overlap sites
polyplot +
geom_point(data = overlap,
aes(x = longitude, y = latitude, shape = dataset),
inherit.aes = FALSE, alpha = 0.6) +
scale_shape_manual(values = c(2, 3, 4))
# ggplot(overlap, aes(x = longitude, y = latitude)) +
# geom_point(aes(shape = dataset, color = Type), alpha = 0.4) +
# scale_shape_manual(values = c(2, 3, 4))
# Combine all aggregated count data and filter for sites in area of overlap
df <- bind_rows(.id = "dataset",
drm = final_counts_ag,
dca = dca_counts_ag,
tt = tt_counts_ag
) %>%
filter(site %in% overlap$site) %>%
select(dataset, site, transect_num, taxon, class, n) %>%
droplevels() %>%
left_join(allsitemd_clean)
write_csv(df, file = "data/processed/all_overlap.csv")2.2. Can ‘Artificial’ be aggregated with ‘Nearshore Ridge Complex’?
‘Artificial’ is only present in DCA and minorly in TT – but bsent from Shedd.
It is in close spatial proximity to Nearshore Ridge Complex – can these be combined?
2.2.1. Test for differences in coral density in DCA survey
# Subset DCA data
dcadf <- dca_counts_ag %>%
left_join(allsitemd_clean) %>%
filter(Type %in% c("Nearshore Ridge Complex", "Inner Reef", "Middle Reef", "Artificial"))
# Fit a Negative Binomial GLM
mod_nb <- MASS::glm.nb(n ~ Type, data = dcadf)
# Generate new data only for existing taxon-size class combinations
newdata_1 <- dcadf %>%
distinct(Type)
# Get predicted values & standard errors (log scale)
preds_nb <- predict(mod_nb, newdata_1, type = "link", se.fit = TRUE)
# Compute both total coral density & taxon-size class-specific densities in one step
results_nb <- newdata_1 %>%
mutate(
fit = exp(preds_nb$fit), # Convert fitted values to response scale
fit_se = exp(preds_nb$fit) * preds_nb$se.fit, # Convert SE using the Delta Method
fit_var = (fit * preds_nb$se.fit)^2, # Variance propagation
fit_lower = exp(preds_nb$fit - 1.96 * preds_nb$se.fit), # Lower CI
fit_upper = exp(preds_nb$fit + 1.96 * preds_nb$se.fit) # Upper CI
)
# Compute total coral density + confidence intervals
total_ci_nb <- results_nb %>%
group_by(Type) %>%
summarize(
total_density = sum(fit),
total_se = sqrt(sum(fit_var)),
lower_95CI = exp(log(total_density) - 1.96 * (total_se / total_density)),
upper_95CI = exp(log(total_density) + 1.96 * (total_se / total_density))
)
knitr::kable(total_ci_nb)| Type | total_density | total_se | lower_95CI | upper_95CI |
|---|---|---|---|---|
| Artificial | 1.7580645 | 0.2269451 | 1.3650635 | 2.264210 |
| Inner Reef | 1.3680352 | 0.1696949 | 1.0727781 | 1.744555 |
| Middle Reef | 0.8573201 | 0.0998316 | 0.6823733 | 1.077120 |
| Nearshore Ridge Complex | 1.9869432 | 0.1764332 | 1.6695542 | 2.364669 |
Highly overlapping confidence intervals for Artifical and Nearshore Ridge Complex, so no difference in coral density.
2.2.2. Test for differences in community composition in DCA survey
library(vegan)
# 1. Create a wide community matrix: rows = site x Type, columns = taxa
comm_matrix <- dcadf %>%
group_by(site, Type, taxon) %>%
summarize(total_n = sum(n), .groups = "drop") %>%
pivot_wider(names_from = taxon, values_from = total_n, values_fill = 0)
# 2. Extract community matrix and metadata
comm_data <- comm_matrix %>% select(-site, -Type)
site_info <- comm_matrix %>% select(site, Type)
# 3. Run NMDS (k = 2 dimensions is standard)
set.seed(123) # for reproducibility
nmds <- metaMDS(comm_data, distance = "bray", k = 2, trymax = 100)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2267843
## Run 1 stress 0.2276397
## Run 2 stress 0.235479
## Run 3 stress 0.2268499
## ... Procrustes: rmse 0.006015082 max resid 0.04502679
## Run 4 stress 0.2325294
## Run 5 stress 0.2437078
## Run 6 stress 0.2276398
## Run 7 stress 0.2356379
## Run 8 stress 0.2315451
## Run 9 stress 0.2297249
## Run 10 stress 0.2314156
## Run 11 stress 0.2374012
## Run 12 stress 0.235405
## Run 13 stress 0.2327499
## Run 14 stress 0.2314157
## Run 15 stress 0.2267804
## ... New best solution
## ... Procrustes: rmse 0.001160623 max resid 0.007358882
## ... Similar to previous best
## Run 16 stress 0.2268231
## ... Procrustes: rmse 0.002309085 max resid 0.018594
## Run 17 stress 0.22764
## Run 18 stress 0.2268497
## ... Procrustes: rmse 0.006092397 max resid 0.04511697
## Run 19 stress 0.2352389
## Run 20 stress 0.2334834
## *** Best solution repeated 1 times# 4. Prepare data for plotting
scores_df <- scores(nmds)$sites %>%
bind_cols(site_info)
# 5. Plot NMDS with ggplot2
ggplot(scores_df, aes(x = NMDS1, y = NMDS2, color = Type)) +
geom_point(size = 3, alpha = 0.8) +
theme_minimal() +
labs(title = "NMDS of Coral Community Structure by Habitat Type",
color = "Habitat Type")
High degree of similarity between Artificial and Nearshore Ridge Complex coral communities.
2.2.3. Combine ‘Artificial’ with ‘Nearshore Ridge Complex’
# Based on these results, combine "Artificial" with "Nearshore Ridge Complex"
df[df$Type == "Artificial", "Type"] <- "Nearshore Ridge Complex"3. Analyze survey data
3.1. Total density
3.1.2. Density per m2 by habitat across surveys
# Compare total corals per square meter across all surveys
# (For DRM, can only include Transect 1 and 2, because 3 and 4 did not record all corals...)
# Add transect size information
df <- df %>%
mutate(transect_num = if_else(is.na(transect_num), 1, transect_num)) %>%
mutate(transect_area_m2 = case_when(
dataset == "drm" ~ 10, # DRM belt transects were 10m
dataset == "tt" ~ 20, # TT belt transects were 20m
dataset == "dca" ~ 30 # DCA belt transects were 30m
)) %>%
mutate(n_per_m2 = n / transect_area_m2)
# Summarize total corals per square meter across all surveys (only T1 and T2 from DRM)
summ <- df %>%
filter(transect_num %in% c(1, 2)) %>%
distinct(dataset, site, transect_num, transect_area_m2) %>% # avoid double-counting transects
group_by(dataset, site) %>%
summarize(
n_transects = n(),
area_per_transect = first(transect_area_m2), # assume it's constant per transect
total_area = n_transects * area_per_transect,
.groups = "drop"
) %>%
left_join(
df %>% group_by(dataset, site, Type) %>% summarize(total_n = sum(n), .groups = "drop"),
by = c("dataset", "site")
) %>%
mutate(total_corals_per_m2 = total_n / total_area)
# Plot by habitat type
ggplot(summ, aes(x = total_corals_per_m2, fill = dataset)) +
facet_grid(Type~.) +
geom_histogram(position = "stack", alpha = 0.5, bins = 30)
3.1.2. Remove lowest density DCA and TT sites
To be extra conservative, filter out tt/dca surveys with total corals per m2 less than the lowest drm survey. In case these would have been considered unsuitable for drm survey, biasing results to higher coral densities in 2024 drm surveys.
summ %>%
select(dataset, site, Type, total_corals_per_m2) %>%
arrange(total_corals_per_m2) %>% print(n = 20)## # A tibble: 241 × 4
## dataset site Type total_corals_per_m2
## <chr> <chr> <chr> <dbl>
## 1 tt Sc14 Nearshore Ridge Complex 0.1
## 2 tt Sc17 Nearshore Ridge Complex 0.1
## 3 tt Sa1 Nearshore Ridge Complex 0.15
## 4 tt Sb6 Nearshore Ridge Complex 0.25
## 5 tt Sc6 Nearshore Ridge Complex 0.25
## 6 dca MI-185-LR Middle Reef 0.267
## 7 dca IN-177-RR Nearshore Ridge Complex 0.3
## 8 tt Sb5 Nearshore Ridge Complex 0.3
## 9 dca IN-194-CP Nearshore Ridge Complex 0.333
## 10 dca MI-36-LR Middle Reef 0.333
## 11 dca MI-37-LR Middle Reef 0.333
## 12 dca MI-181-LR Middle Reef 0.367
## 13 dca MI-39-LR Middle Reef 0.4
## 14 drm 3072 Nearshore Ridge Complex 0.45
## 15 tt Nf2 Middle Reef 0.45
## 16 tt Sb3 Nearshore Ridge Complex 0.45
## 17 tt Sc10 Nearshore Ridge Complex 0.45
## 18 tt Sc12 Nearshore Ridge Complex 0.45
## 19 dca IN-2-AR Nearshore Ridge Complex 0.467
## 20 dca IN-87-RR Nearshore Ridge Complex 0.5
## # ℹ 221 more rowslowsites <- summ %>% filter(total_corals_per_m2 < 0.45) %>% pull(site)
dff <- df %>%
filter(!site %in% lowsites)3.1.3. Remove lowest abundance coral taxa
These will be problematic for statistical models
# Drop taxa with very few observations
sppcounts <- dff %>%
group_by(taxon) %>%
summarize(total = sum(n), .groups = "drop") %>%
arrange(total)
sppcounts## # A tibble: 23 × 2
## taxon total
## <chr> <int>
## 1 HCUC 0
## 2 MANG 0
## 3 PAME 0
## 4 FFRA 1
## 5 ODIF 1
## 6 OROB 1
## 7 SCOL 1
## 8 ORBI 5
## 9 MUSS 12
## 10 MYCE 19
## # ℹ 13 more rowsdff <- dff %>%
filter(taxon %in% filter(sppcounts, total >= 5)$taxon)Taxa with less than 5 total observations were filtered out, which included: HCUC, MANG, PAME, FFRA, ODIF, OROB, SCOL
3.2. Taxon/size class density
3.2.1. By survey (overall)
Negative binomial model with taxon-sizeclass
dff <- dff %>% mutate(taxclass = interaction(taxon, class, sep = ""))
# Fit a Negative Binomial GLM
mod_nb <- MASS::glm.nb(n ~ dataset * taxclass + offset(log(transect_area_m2)), data = dff)
# Generate new data only for existing taxon-size class combinations
newdata_1 <- dff %>%
distinct(dataset, taxon, class, taxclass) %>% # Keep only observed taxon-class combinations
mutate(transect_area_m2 = 1) # Set area to 1 for density predictions on a per m2 basis
# Get predicted values & standard errors (log scale)
preds_nb <- predict(mod_nb, newdata_1, type = "link", se.fit = TRUE)
# Compute both total coral density & taxon-size class-specific densities in one step
results_nb <- newdata_1 %>%
mutate(
fit = exp(preds_nb$fit), # Convert fitted values to response scale
fit_se = exp(preds_nb$fit) * preds_nb$se.fit, # Convert SE using the Delta Method
fit_var = (fit * preds_nb$se.fit)^2, # Variance propagation
fit_lower = exp(preds_nb$fit - 1.96 * preds_nb$se.fit), # Lower CI
fit_upper = exp(preds_nb$fit + 1.96 * preds_nb$se.fit) # Upper CI
)
# Compute total coral density + confidence intervals
total_ci_nb <- results_nb %>%
group_by(dataset) %>%
summarize(
total_density = sum(fit),
total_se = sqrt(sum(fit_var)),
lower_95CI = exp(log(total_density) - 1.96 * (total_se / total_density)),
upper_95CI = exp(log(total_density) + 1.96 * (total_se / total_density))
)
knitr::kable(total_ci_nb)| dataset | total_density | total_se | lower_95CI | upper_95CI |
|---|---|---|---|---|
| dca | 1.692557 | 0.1035547 | 1.501287 | 1.908195 |
| drm | 2.657382 | 0.1726744 | 2.339606 | 3.018321 |
| tt | 1.693507 | 0.1011996 | 1.506331 | 1.903941 |
# 2.66 total corals per m2 in drm (2.63 if you filter out the further south sites using 'overlap_lat')
# 1.69 in DCA
# 1.69 in TT
# Extract and plot taxon-size-class-specific densities
taxclass_ci_nb <- results_nb %>%
dplyr::select(dataset, taxon, class, taxclass, fit, fit_se, fit_lower, fit_upper)
# Compute total density per taxon (summing over size classes)
fitted_taxon_nb <- taxclass_ci_nb %>%
group_by(dataset, taxon) %>%
summarize(
fit = sum(fit), # Sum densities across size classes
fit_se = sqrt(sum(fit_se^2)), # Correct SE propagation (variance summation)
log_fit = log(fit), # Log-transform fit for proper CI computation
log_se = fit_se / fit, # Approximate log-scale standard error
fit_lower = exp(log_fit - 1.96 * log_se), # Compute lower CI on log scale
fit_upper = exp(log_fit + 1.96 * log_se) # Compute upper CI on log scale
) %>%
ungroup() %>%
mutate(taxon = fct_reorder(taxon, fit), class = "Total") # Reorder taxa by abundance
# fitted_combined_nb <- bind_rows(taxclass_ci_nb, fitted_taxon_nb) %>%
# mutate(taxon = factor(taxon, levels = levels(fitted_taxon_nb$taxon)))
# Combine taxclass and taxon totals, and remove taxa/classes not observed in a dataset
fitted_combined_nb <- taxclass_ci_nb %>% right_join(dff %>% filter(n > 0) %>% distinct(dataset, taxon, class)) %>%
bind_rows(fitted_taxon_nb %>% right_join(dff %>% filter(n > 0) %>% distinct(dataset, taxon))) %>%
mutate(taxon = factor(taxon, levels = levels(fitted_taxon_nb$taxon)))
# Plot size classes and taxon totals
ggplot(fitted_combined_nb, aes(x = taxon, y = fit, color = dataset, shape = dataset)) +
geom_point(aes(size = class),
position = position_dodge(width = 0.2), alpha = 0.6) +
geom_errorbar(aes(ymin = fit_lower, ymax = fit_upper),
width = 0, position = position_dodge(width = 0.2), alpha = 0.6) +
scale_y_log10(limits = c(1e-4, 5)) +
scale_size_manual(values = c("Total" = 3, ">4cm" = 2, "<4cm" = 1)) + # Larger points for totals
#scale_shape_manual(values = c("Total" = 15, ">4cm" = 16, "<4cm" = 16)) + # Different shape for totals
coord_flip() +
labs(y = "Estimated Coral Density (per m²)", x = "Taxon",
color = "Dataset", shape = "Size class", size = "Size class") +
theme_minimal() +
labs(title = "NB Model")
# Plot just taxon totals
fitted_combined_nb %>%
filter(class == "Total") %>%
ggplot(aes(x = taxon, y = fit, color = dataset)) + #, shape = dataset)) +
geom_point(position = position_dodge(width = 0.2), alpha = 0.6) +
geom_errorbar(aes(ymin = fit_lower, ymax = fit_upper),
width = 0, position = position_dodge(width = 0.2), alpha = 0.6) +
scale_y_log10(limits = c(1e-4, 5)) +
#scale_shape_manual(values = c("Total" = 15, ">4cm" = 16, "<4cm" = 16)) + # Different shape for totals
coord_flip() +
labs(y = "Estimated Coral Density (per m²)", x = "Taxon",
color = "Dataset", shape = "Size class", size = "Size class") +
theme_minimal() +
labs(title = "NB Model")
3.2.2. By habitat by survey
Negative binomial model with taxon and habitat Type
dfftaxon <- dff %>%
group_by(dataset, Type, site, transect_num, transect_area_m2, taxon) %>%
summarize(n = sum(n)) %>%
ungroup()
# Fit a Negative Binomial GLM
# Note: also tried pscl::zeroinfl and was nearly identical to glm.nb results
mod_nb <- MASS::glm.nb(n ~ dataset * Type * taxon + offset(log(transect_area_m2)), data = dfftaxon)
# Generate new data only for existing taxon-size class combinations
newdata_1 <- dfftaxon %>%
distinct(dataset, Type, taxon) %>%
mutate(transect_area_m2 = 1) # Set area to 1 for density predictions on a per m2 basis
# Get predicted values & standard errors (log scale)
preds_nb <- predict(mod_nb, newdata_1, type = "link", se.fit = TRUE)
newdata_1 <- bind_cols(newdata_1, preds_nb) %>% mutate(fit = exp(fit)) %>%
mutate(taxon = fct_reorder(taxon, fit), class = "Total") # Reorder taxa by abundance
# Compute both total coral density & taxon-size class-specific densities in one step
results_nb <- newdata_1 %>%
mutate(
fit = exp(preds_nb$fit), # Convert fitted values to response scale
fit_se = exp(preds_nb$fit) * preds_nb$se.fit, # Convert SE using the Delta Method
fit_var = (fit * preds_nb$se.fit)^2, # Variance propagation
fit_lower = exp(preds_nb$fit - 1.96 * preds_nb$se.fit), # Lower CI
fit_upper = exp(preds_nb$fit + 1.96 * preds_nb$se.fit) # Upper CI
)
# Compute total coral density + confidence intervals
total_ci_nb <- results_nb %>%
group_by(dataset, Type) %>%
summarize(
total_density = sum(fit),
total_se = sqrt(sum(fit_var)),
lower_95CI = exp(log(total_density) - 1.96 * (total_se / total_density)),
upper_95CI = exp(log(total_density) + 1.96 * (total_se / total_density))
)
knitr::kable(total_ci_nb)| dataset | Type | total_density | total_se | lower_95CI | upper_95CI |
|---|---|---|---|---|---|
| dca | Inner Reef | 1.413636 | 0.1579489 | 1.1356100 | 1.759731 |
| dca | Middle Reef | 1.015873 | 0.0980266 | 0.8408167 | 1.227376 |
| dca | Nearshore Ridge Complex | 2.031667 | 0.1768596 | 1.7129817 | 2.409640 |
| drm | Inner Reef | 2.345674 | 0.2349990 | 1.9274773 | 2.854605 |
| drm | Middle Reef | 2.116667 | 0.3498921 | 1.5308879 | 2.926588 |
| drm | Nearshore Ridge Complex | 2.872303 | 0.2602578 | 2.4049269 | 3.430509 |
| tt | Inner Reef | 1.961539 | 0.2889501 | 1.4696226 | 2.618110 |
| tt | Middle Reef | 1.531250 | 0.1525842 | 1.2595779 | 1.861518 |
| tt | Nearshore Ridge Complex | 1.703750 | 0.1561501 | 1.4236103 | 2.039016 |
# Extract and plot taxon-size-class-specific densities
taxclass_ci_nb <- results_nb %>%
dplyr::select(dataset, Type, taxon, fit, fit_se, fit_lower, fit_upper)
# Compute total density per taxon (summing over size classes)
fitted_taxon_nb <- taxclass_ci_nb %>%
group_by(dataset, Type, taxon) %>%
summarize(
fit = sum(fit), # Sum densities across size classes
fit_se = sqrt(sum(fit_se^2)), # Correct SE propagation (variance summation)
log_fit = log(fit), # Log-transform fit for proper CI computation
log_se = fit_se / fit, # Approximate log-scale standard error
fit_lower = exp(log_fit - 1.96 * log_se), # Compute lower CI on log scale
fit_upper = exp(log_fit + 1.96 * log_se) # Compute upper CI on log scale
) %>%
ungroup() %>%
mutate(taxon = fct_reorder(taxon, fit), class = "Total") # Reorder taxa by abundance
fitted_combined_nb <- bind_rows(taxclass_ci_nb, fitted_taxon_nb) %>%
mutate(taxon = factor(taxon, levels = levels(fitted_taxon_nb$taxon)))
# Combine taxclass and taxon totals, and remove taxa/classes not observed in a dataset
fitted_combined_nb <- taxclass_ci_nb %>% right_join(dff %>% filter(n > 0) %>% distinct(dataset, taxon, Type)) %>%
bind_rows(fitted_taxon_nb %>% right_join(dff %>% filter(n > 0) %>% distinct(dataset, taxon, Type))) %>%
mutate(taxon = factor(taxon, levels = levels(fitted_taxon_nb$taxon)))
ggplot(fitted_combined_nb, aes(x = taxon, y = fit, color = dataset)) +
geom_point(aes(size = class),
position = position_dodge(width = 0.2), alpha = 0.6) +
geom_errorbar(aes(ymin = fit_lower, ymax = fit_upper),
width = 0, position = position_dodge(width = 0.2), alpha = 0.6) +
facet_wrap(~Type) +
scale_y_log10(limits = c(1e-4, 5)) +
#scale_size_manual(values = c("Total" = 4, ">4cm" = 2.5, "<4cm" = 1)) + # Larger points for totals
#scale_shape_manual(values = c("Total" = 15, ">4cm" = 16, "<4cm" = 16)) + # Different shape for totals
coord_flip() +
labs(y = "Estimated Coral Density (per m²)", x = "Taxon",
color = "Dataset", shape = "Size class", size = "Size class") +
theme_minimal() +
labs(title = "NB Model")
4. Total corals in impact zones
4.1. Calculate area of each habitat type in each impact zone
impact_zones <- st_read("data/Impact_zones.kml") %>%
st_zm(drop = TRUE, what = "ZM") %>%
st_make_valid() %>%
st_transform(st_crs(polygons_clean))## Reading layer `Impact_zones' from data source
## `/Users/rosscunning/Projects/PENIP/data/Impact_zones.kml' using driver `KML'
## Simple feature collection with 9 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: -80.10538 ymin: 26.08289 xmax: -80.08298 ymax: 26.1069
## z_range: zmin: -21.40724 zmax: 6.930363e-14
## Geodetic CRS: WGS 84impact_zones <- impact_zones %>%
rename(ImpactZone = Name) # or whatever column contains zone names
library(ggplot2)
impact_zones_plot <- impact_zones %>%
mutate(linetype_group = "Impact Zone")
ggplot() +
# Habitat polygons
geom_sf(data = polygons_clean, aes(fill = Type), color = "black", size = 0.2, alpha = 0.6) +
# Impact zones with dummy linetype for legend
geom_sf(data = impact_zones_plot, aes(linetype = linetype_group),
fill = NA, color = "black", linewidth = 0.6, show.legend = TRUE) +
# Color scales
scale_fill_brewer(palette = "Set3", na.value = "gray80") +
scale_linetype_manual(name = "", values = c("Impact Zone" = "dashed")) +
# Theme and layout
theme_minimal() +
labs(title = "Habitat Polygons within Impact Zones", fill = "Type") +
theme(legend.position = "right") +
xlim(-80.11, -80.079) +
ylim(26.08, 26.11)
# Spatial intersection of habitat polygons with impact zones
habitat_in_zones <- st_intersection(polygons_clean, impact_zones)
# Use a projected CRS for accurate area (e.g., UTM Zone 17N for South Florida)
habitat_in_zones_proj <- habitat_in_zones %>%
st_transform(32617) %>%
mutate(area_m2 = st_area(geometry))
# Adjust column names depending on your Impact Zones KML
area_summary <- habitat_in_zones_proj %>%
st_drop_geometry() %>%
group_by(Type, ImpactZone) %>%
summarize(total_area_m2 = sum(as.numeric(area_m2)), .groups = "drop")
area_summary2 <- area_summary %>%
filter(Type %in% c("Nearshore Ridge Complex", "Artificial", "Inner Reef", "Middle Reef"))4.2. Calculate corals in each impact zone
area_summary2## # A tibble: 27 × 3
## Type ImpactZone total_area_m2
## <chr> <chr> <dbl>
## 1 Artificial Scenario 2, 5.1-10 cm 11862.
## 2 Artificial Scenario 2, > 10 cm 3404.
## 3 Artificial Scenario 4, 0.1-0.5 cm 2310.
## 4 Artificial Scenario 4, 0.51-1 cm 60506.
## 5 Artificial Scenario 4, 1.1-5 cm 288364.
## 6 Artificial Side Slopes 3847.
## 7 Inner Reef Channel 517.
## 8 Inner Reef Scenario 2, 5.1-10 cm 11048.
## 9 Inner Reef Scenario 2, > 10 cm 12133.
## 10 Inner Reef Scenario 4, 0.1-0.5 cm 415915.
## # ℹ 17 more rowsfitted_combined_nb## # A tibble: 238 × 10
## dataset Type taxon fit fit_se fit_lower fit_upper log_fit log_se class
## <chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 dca Inner… AGAR 0.0106 0.00445 0.00466 0.0242 NA NA <NA>
## 2 dca Inner… DSTO 0.0182 0.00621 0.00931 0.0355 NA NA <NA>
## 3 dca Inner… EFAS 0.00303 0.00221 0.000724 0.0127 NA NA <NA>
## 4 dca Inner… FAVI 0.0288 0.00844 0.0162 0.0512 NA NA <NA>
## 5 dca Inner… MADR 0.00152 0.00154 0.000207 0.0111 NA NA <NA>
## 6 dca Inner… MCAV 0.192 0.0391 0.129 0.287 NA NA <NA>
## 7 dca Inner… MEAN 0.0273 0.00813 0.0152 0.0489 NA NA <NA>
## 8 dca Inner… MMEA 0.00152 0.00154 0.000207 0.0111 NA NA <NA>
## 9 dca Inner… MYCE 0.00152 0.00154 0.000207 0.0111 NA NA <NA>
## 10 dca Inner… PORI 0.100 0.0220 0.0649 0.154 NA NA <NA>
## # ℹ 228 more rowstotals <- full_join(fitted_combined_nb, area_summary2, by = "Type") %>%
mutate(tot_lower = fit_lower * total_area_m2,
tot_upper = fit_upper * total_area_m2) %>%
select(dataset, Type, taxon, ImpactZone, tot_lower, tot_upper)
res <- totals %>%
drop_na(dataset) %>%
group_by(dataset, ImpactZone) %>%
summarize(totlow = sum(tot_lower),
totup = sum(tot_upper)) %>%
pivot_wider(names_from = dataset, values_from = c(totlow, totup))
res2 <- res %>%
# Add totals
add_row(ImpactZone = "Total",
!!!res %>% summarise(across(where(is.numeric), sum, na.rm = TRUE)) %>% as.list()) %>%
# Format ranges
mutate(across(where(is.numeric), ~ sprintf("%.1e", .))) %>%
unite(dca, totlow_dca, totup_dca, sep = " — ") %>%
unite(tt, totlow_tt, totup_tt, sep = " — ") %>%
unite(drm, totlow_drm, totup_drm, sep = " — ")
knitr::kable(res2, caption = "Total number of corals in each Impact Zone for the Nearshore Ridge Complex, Inner, and Middle Reefs. Note: DOES NOT INCLUDE Outer Reef.")| ImpactZone | dca | drm | tt |
|---|---|---|---|
| Channel | 4.2e+04 — 1.0e+05 | 7.7e+04 — 2.4e+05 | 6.3e+04 — 1.4e+05 |
| Scenario 2, 5.1-10 cm | 1.4e+05 — 2.5e+05 | 2.1e+05 — 4.0e+05 | 1.3e+05 — 2.7e+05 |
| Scenario 2, > 10 cm | 1.7e+05 — 3.0e+05 | 2.4e+05 — 4.7e+05 | 1.5e+05 — 3.2e+05 |
| Scenario 4, 0.1-0.5 cm | 4.5e+06 — 7.9e+06 | 6.5e+06 — 1.2e+07 | 4.0e+06 — 8.6e+06 |
| Scenario 4, 0.51-1 cm | 6.3e+05 — 1.1e+06 | 9.3e+05 — 1.8e+06 | 5.8e+05 — 1.3e+06 |
| Scenario 4, 1.1-5 cm | 6.7e+05 — 1.3e+06 | 1.0e+06 — 2.1e+06 | 6.6e+05 — 1.5e+06 |
| Side Slopes | 5.6e+04 — 9.6e+04 | 8.1e+04 — 1.6e+05 | 4.9e+04 — 9.8e+04 |
| Total | 6.2e+06 — 1.1e+07 | 9.1e+06 — 1.7e+07 | 5.6e+06 — 1.2e+07 |